#include "mcmcSymphony.h"
//#include "IParsingDoc.h"
//#include "Resource.h"
//#include "MainFrm.h"
#include "ddmcmc.h"

//extern CIParsingDoc *pDocActive;
extern DDMCMC *DDMCMCActive;

#define DN_SWITCH		0
#define DN_DIFFUSION	1

extern long long tFunction;
long long tStart;

/////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// operations
MCMC_Conductor::MCMC_Conductor(MCMC_Map_2Layers &map, mcmcMatrix<int> partition_maps1[3], int num[],Raster_Color<uchar> &raster,
							   mcmcMatrix<int> &mx_init_clst, const MCMC_SegSetting& setting, char *pfile_name)
{
	// initialization
	int	num_region_segmented[3];
	rf_dynamics.Set(0,NUM_DYNAMIC-1,1);	// split or merge
	rf_dynamics.InitValue(1.0);
	rf_dynamics.SetCdf();

	bcolor = setting.bcolor;

	pmap = &map;

	respFilter.Set(raster,IsColorImage()); // Time Consuming: 100ms
	respFilter.SetAllFilterBanks(); // Time Consuming: 500ms
	

	pmap->SetFilterResp(respFilter);
	step = -1;

	// create particles
	pmcmc_particle_gauss = NULL;


	
	if (!IsColorImage())
	{
		pmcmc_particle_gauss = new MCMC_Particle_Intensity();
		pmcmc_particle_gauss->Set(respFilter,mx_init_clst,pfile_name);
	}
	else
	{
#ifdef LUV_COLOR_MODEL
		pmcmc_particle_gauss = new MCMC_Particle_3DColor2();
#else
		pmcmc_particle_gauss = new MCMC_Particle_3DColor();
#endif
		pmcmc_particle_gauss->Set(respFilter,mx_init_clst,pfile_name); // time consuming: 2100ms
	}
	

	// use the clustring map as edge map proposal at scale 1
/*
	Raster<int> ras_cluster;
	ras_cluster = pmcmc_particle_gauss->mx_map_temp;
	num[0] = ras_cluster.ResetLabels(false);
	for (int j=0; j<raster.rows(); j++)
		for (int i=0; i<raster.cols(); i++)
			seg_regions[0](j,i) = ras_cluster(j,i);
*/
	
	// initialization of diffusion class
	mcmc_diffusion.Set(*pmap,raster,pmcmc_particle_gauss);
	model_supervisor.Set(*pmap,*pmcmc_particle_gauss,respFilter,setting);

	for (int i=0; i<3; i++)
	{
		partition_maps[i] = &partition_maps1[i];
		num_region_segmented[i] = num[i];
		// cout << "num: " << num[i] << endl; // Siyuan
	}
	
	// preparation
	model_supervisor.InitRegionModels();
	
	// compute the energy
	pmap->ComputeTotalEnergy(RECOMPUTE_ENERGY_SHAPE);

	// map buffer
	pmap_buf = new MCMC_Map_2Layers();
	assert(pmap_buf != NULL);
	
	// constructions of dynamics
	dynamics[0] = new MCMC_Dn_Switch(*pmap, respFilter, model_supervisor,step);
	dynamics[1] = new MCMC_Dn_Diffusion(*pmap, respFilter, mcmc_diffusion,model_supervisor,step);
	sw_coordinator.Set(*map.LayerRegion(), model_supervisor);
	pdn_change_underling = new MCMC_Dn_RepartitionARegions(*pmap, respFilter, model_supervisor, partition_maps1);
	
}

MCMC_Conductor::~MCMC_Conductor()
{
	pmap = NULL;

	// free objects
	for (int i=0; i<NUM_DYNAMIC; i++)
	{
		if (dynamics[i] != NULL)
			delete dynamics[i];
		dynamics[i] = NULL;
	}

	if (!IsColorImage())
	{
		if (pmcmc_particle_gauss != NULL)
		{
			MCMC_Particle_Intensity *par_temp;
			par_temp = (MCMC_Particle_Intensity *)pmcmc_particle_gauss;
			delete par_temp;
			pmcmc_particle_gauss = NULL;
		}
	}
	else
	{
		if (pmcmc_particle_gauss != NULL)
		{
#ifdef 	LUV_COLOR_MODEL
			MCMC_Particle_3DColor2 *par_temp;
			par_temp = (MCMC_Particle_3DColor2 *)pmcmc_particle_gauss;
#else
			MCMC_Particle_3DColor *par_temp;
			par_temp = (MCMC_Particle_3DColor *)pmcmc_particle_gauss;
#endif
			delete par_temp;
			pmcmc_particle_gauss = NULL;
		}
	}

	// map buffer
	if (pmap_buf != NULL)
		delete pmap_buf;
	pmap_buf = NULL;

	if (pdn_change_underling != NULL)
		delete pdn_change_underling;
	pdn_change_underling = NULL;
}

void MCMC_Conductor::SetModelPointers(void)
{
	if (LayerRegion()!=NULL)
	{
		int			id_region;
		MCMC_AtomicRegion *pregion;
		MCMC_Model	*pmodel;

		id_region = LayerRegion()->GetFirstAtomicRegion();
		while (id_region>=0)
		{
			pregion = LayerRegion()->GetAtomicRegion(id_region);
			assert(pregion!=NULL);
			pmodel = pregion->Model();

			if (pmodel!=NULL)
			{
				pmodel->SetImage(respFilter);
				pmodel->SetParticel(*model_supervisor.Particle(pmodel->Type()));
			}
			id_region = LayerRegion()->GetNextAtomicRegion();
		}
	}
}

double MCMC_Conductor::UpdateEnergy(const int step_energy)
{
	double energy_temp;

	pmap->ComputeTotalEnergy(DONOT_RECOMPUTE_ENERGY_SHAPE);
	energy_temp = pmap->EnergyTotal();

	if (step-1>=0 && vc_energy_total(step-1)>1e2 && energy_temp-vc_energy_total(step-1)>7e3)
		energy_temp = vc_energy_total(step-1)+(energy_temp-vc_energy_total(step-1))/100.0;

	vc_energy_total(step) = energy_temp;

	return energy_temp;
}

void MCMC_Conductor::SaveEnergy(char *file_name)
{
	FILE *file;
	file = fopen(file_name, "w+");
	if (file != NULL)
	{
		for (int i=0; i<vc_energy_total.size(); i++)
			fprintf(file,"%.3lf\n", vc_energy_total(i));
		fclose(file);
	}
}

void MCMC_Conductor::ChangeAllGraphPointers(MCMC_Map_2Layers &map)
{
	pmap = &map;

	mcmc_diffusion.ChangeMapPointer(*pmap);
	for(int i=0; i<NUM_DYNAMIC; i++)
	{
		if (dynamics[i]!=NULL)
			dynamics[i]->ChangeMapPointer(*pmap);
	}
	model_supervisor.ChangeMapPointer(*pmap);
	scale_factor = LayerRegion()->scale_factor;
}

void MCMC_Conductor::SetScaleFactor(const double sc)
{
	LayerRegion()->scale_factor = sc;
	scale_factor = LayerRegion()->scale_factor;
}


bool MCMC_Conductor::OneStep(void)
{
	bool bAccept=false;

	MCMC_AtomicRegion *pregion;
	int id_region = LayerRegion()->GetFirstAtomicRegion();
	while (id_region>=0)
	{
		pregion = LayerRegion()->GetAtomicRegion(id_region);
		id_region = LayerRegion()->GetNextAtomicRegion();
	}

	// select a dynamic to go
	dn = (int)rf_dynamics.Sampling();
//ztu_temp
//dn = 0; // Siyuan!

	// make a proposal
	dynamics[dn]->Propose();
	// accept?
	bAccept = dynamics[dn]->TryAccept();
	// do the corresponding action
	dynamics[dn]->DoAction();
	//pDocActive->DisplayInfo(dynamics[dn]->Info); Jiadi Yang

	if (bAccept)
		LayerRegion()->ComputeTotalEnergy(DONOT_RECOMPUTE_ENERGY_SHAPE);
/*
	if (pmap->HasCurveLayer())
	{
		if (pmap->GetT() < 1500.0)
		{
			if (pmap->Step()%10==0)
				pcoordinator_curve->OneStep(true);
			else
				pcoordinator_curve->OneStep(false);

			for (int k=0; k<1; k++)
				if (pcoordinator_curve->Dn()!=6)
				  pcoordinator_curve->OneStep(false);

			if (pmap->GetT() < 100.0 && pcoordinator_curve->Dn()!=6)
				pcoordinator_curve->OneStep(false);
			if (pmap->GetT() < 10.0 && pcoordinator_curve->Dn()!=6)
				pcoordinator_curve->OneStep(false);
		}
	}
*/
	// update current energy
	UpdateEnergy();
	// update the energy curve in the doc

	// Siyuan
	Mat mSyn, mBoundaries;
	if (step % 20 == 0){
		mSyn = DDMCMCActive->synthesize();
		//mBoundaries = DDMCMCActive->boundaries();
	}
	imshow("Synthesized in process (OneStep)", mSyn);
	//imshow("Boundries in process (OneStep)", mBoundaries);
	waitKey(1);


	/* Jiadi Yang
	pDocActive->SetStep(step, vc_energy_total(step)); 
	if (step%20==0)
		pDocActive->ShowSynthesizedImage();
 	pDocActive->RepaintEnergyViews();
	*/
	return bAccept;
}


void MCMC_Conductor::Go(char *fname, const MCMC_SegSetting *psetting)
{
	int					total;
	double				e_temp=1e10;
	char				str[200];
	//CIParsingDoc			*pDoc; Jiadi Yang

	// to the number of steps
	if (psetting==NULL)
		total = 2000;
	else
		total = psetting->nstep_to;

	// from the number of steps
	if (psetting==NULL)
		step = 0;
	else
		step = psetting->nstep_from;

	//step = 5001;
	total = 5000;

FILE	*pfile;
pfile = fopen("rand.txt","w+t");
fclose(pfile);

	vc_energy_total.SetDimension(total,1);
	vc_energy_total.InitValue(0.0);

	// pDocActive->SetEnergy(vc_energy_total, step);  Jiadi Yang

	SetModelPointers();

	// set the scale factor
	if (psetting==NULL)
		SetScaleFactor(7.0);
	else
		SetScaleFactor(psetting->dscale);

	if (psetting!= NULL && psetting->bload_data)
	{
		pmap->ReadFromFile(psetting->pfile_name);
		SetScaleFactor(pmap->LayerRegion()->scale_factor);
	}

	LayerRegion()->CreateTopologyRecords();
	model_supervisor.FitAllRegionParms();

	// compute the total energy
	pmap->ComputeTotalEnergy(RECOMPUTE_ENERGY_SHAPE);
	pmap->Reconstruct();

	//pDocActive->ShowSynthesizedImage(); Jiadi Yang

	while (step<total)
	{
		pmap->SetStep(step);

		// one dynamic
		pfile = fopen("rand.txt","at");
		fprintf(pfile, "iteration=%d, rand seed=%u\n", step, ztu_seed());
		fclose(pfile);

		OneStep();

		step++;
	}

	// save results
	Raster<int>		image;
	Raster<uchar>	raster_boundaries;
	int				ii,jj;

	raster_boundaries.SetDimension(respFilter.rows(),respFilter.cols());
	image.SetDimension(respFilter.rows(),respFilter.cols());
	image.InitValue(0);

	LayerRegion()->DrawRegions(image);
	// show the boundaries only
	for (jj=0; jj<image.rows(); jj++)
		for (ii=0; ii<image.cols(); ii++)
		{
			if (image(jj,ii)>0)
				raster_boundaries(jj,ii) = 0;
			else
				raster_boundaries(jj,ii) = 255;
		}

	sprintf(str,"boundaries_%.1f", LayerRegion()->scale_factor);
	//pDoc = ((CMainFrame*)AfxGetMainWnd())->NewImageDoc(raster_boundaries,10,10,str); // Jiadi Yang, commented but not yet rewritten
}

void MCMC_Conductor::OneStepNew(void)
{
	int		ndynamics;
	
	//the kernel for SW-Cuts
	rf_dynamics_new.Set(-1, NUM_DYNAMIC-1, 1);
	rf_dynamics_new.InitValue(0.0);

	rf_dynamics_new.SetOneBin(-1, 0.7);
	rf_dynamics_new.SetOneBin(DN_SWITCH, 0.3);
	if (pmap->Step() >= 500)
		rf_dynamics_new.SetOneBin(DN_DIFFUSION, 0.002);
	rf_dynamics_new.normalize();

	ndynamics = rf_dynamics_new.SamplingIntNonSetCdf();

	if (ndynamics == -1)
	{
		//SW cuts
		sw_coordinator.OneStep();
		
		//pDocActive->DisplayInfo(sw_coordinator.Info); Jiadi Yang
			
		sw_coordinator.ClearShowProposal();
	}
	else
	{
		/*
		if (ndynamics == 1){
			int fooDebug = 0;
			cout << "dynamic: diffusion" << endl;
		}
		*/
		// make a proposal
		dynamics[ndynamics]->Propose();
		// accept?
		dynamics[ndynamics]->TryAccept();
		// do the corresponding action
		dynamics[ndynamics]->DoAction();

		//pDocActive->DisplayInfo(sw_coordinator.Info); Jiadi Yang
	}
}

bool CheckCheck(MCMC_Layer_Regions *player_region)
{
	// to synthesize
	LinkList<int>		list_aregions,list_indexid;
	MCMC_AtomicRegion	*paregion;
	bool				bgo,bgo2;
	int indexid,i,j;
	MCMC_RegionElement	*pregion;
	
	mcmcMatrix<uchar> mx_tt;

	mx_tt.SetDimension(player_region->height, player_region->width);
	mx_tt.InitValue(0);

	player_region->GetAllTopNodeRegions(list_indexid);
	bgo = list_indexid.GetFirst(indexid);
	while (bgo)
	{
		pregion = player_region->GetRegionFromIndexId(indexid);
			
			paregion = player_region->GetAtomicRegion(pregion->GetFirstAtomicRegion());
			while (paregion != NULL)
			{
				bgo2 = paregion->GetFirstPixel(j, i);
				while (bgo2)
				{
					
					mx_tt(j ,i)=1;
					bgo2 = paregion->GetNextPixel(j, i);

				}
				paregion = player_region->GetAtomicRegion(pregion->GetNextAtomicRegion());
			}
			bgo = list_indexid.GetNext(indexid);
		}
	for (j=0; j<mx_tt.rows(); j++)
		for (i=0; i<mx_tt.cols(); i++)
		{
			if (mx_tt(j, i)==0)
				return false;
		}
	return true;
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
void MCMC_Conductor::GoNew(char *fname, const MCMC_SegSetting &setting)
{
	int		total;

	step = setting.nstep_from;
	total = setting.nstep_to;

	vc_energy_total.SetDimension(total,1);
	vc_energy_total.InitValue(0.0);
 
	 //pDocActive->SetEnergy(vc_energy_total,step); Jiadi Yang

//pmap->ReadFromFile("a.dat");


	// initialization
	SetModelPointers();
	// set the scale factor
	SetScaleFactor(setting.dscale);
	LayerRegion()->CreateTopologyRecords(); // Time consuming: 1300ms

	model_supervisor.FitAllRegionParms();

	// initialize the affinity map
	LayerRegion()->CreateAffinityMap(respFilter);

	char str_data_name[200];

	sprintf(str_data_name, "%s.data", fname);

	// compute the total energy
	pmap->ComputeTotalEnergy(RECOMPUTE_ENERGY_SHAPE);
	//pDocActive->ShowSynthesizedImage(); Jiadi Yang

	// Siyuan Qi
	//Mat mSyn = DDMCMCActive->synthesize();//, mBoundaries = DDMCMCActive->boundaries();
	//imshow("Synthesized in process", mSyn);
	//imshow("Boundries in process", mBoundaries);
	//waitKey(1);

	while (step<total)
	{
		// Siyuan
		//if (step % 20 == 0){
			//mSyn = DDMCMCActive->synthesize();
			//mBoundaries = DDMCMCActive->boundaries();
		//}
		//imshow("Synthesized in process", mSyn);
		//imshow("Boundries in process", mBoundaries);
		//waitKey(1);
		// Jiadi
		if (setting.visualize){
			Raster<int>			mx_boundary;
			Mat mBoundary = Mat(respFilter.rows(), respFilter.cols(), CV_8U);
			mx_boundary.SetDimension(respFilter.rows(), respFilter.cols());
			mx_boundary.InitValue(0);
			pmap->LayerRegion()->DrawRegions(mx_boundary);
			for (int j = 0; j < mBoundary.rows; j++){
				for (int i = 0; i<mBoundary.cols; i++){
					if (mx_boundary(j, i)>0){
						mBoundary.at<uchar>(j, i) = 255;
					}
					else {
						mBoundary.at<uchar>(j, i) = 0;
					}
				}
			}
			imshow("Boundaries in Process", mBoundary);
			waitKey(1);
		}
		// End Jiadi

		//cout << step << endl;
		pmap->SetStep(step);

		OneStepNew();

		UpdateEnergy();
		// update the energy curve in the doc
		//pDocActive->SetStep(step,vc_energy_total(step)); Jiadi Yang

		if (step==2000)
			pdn_change_underling->Propose();
		
		//if (step%50 == 0)
			//pDocActive->OnRepaintAllViews(); Jiadi Yang

		step++;
	}

	Raster<int>		ras_labels;
	Raster<uchar>	ras_aa;

	sprintf(str_data_name, "%s_labels.tif", fname);
	LayerRegion()->LabelMapTopRegions(ras_labels);
	CopyRaster(ras_aa, ras_labels);
	SaveRaster2ImageFile(ras_aa, str_data_name);
}

